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Simplified Solar Modulation Model of Inner Trapped Belt Proton Flux as a 

Function of Atmospheric Density 


ABSTRACT 


No simple algorithm seems to exist for calculating proton fluxes and lifetimes in 
the Earth’s inner, trapped radiation belt throughout the solar cycle. Most models of the 
inner trapped belt in use depend upon APS which only describes the radiation 
environment at solar maximum and solar minimum in Cycle 20. One exception is 
NOAAPRO which incorporates flight data from the TIROS/NOAA polar orbiting 
spacecraft. The present study discloses yet another, simple formulation for 
approximating proton fluxes at any time in a given solar cycle, in particular between solar 
maximum and solar minimum. It is derived from AP8 using a regression algorithm 
technique from nuclear physics. From flux and its time integral fluence, one can then 
approximate dose rate and its time integral dose. It has already been published in this 
journal that the absorbed dose rate, D, in the trapped belts exhibits a power law 
relationship, D - Ap' n , where A is a constant, p is the atmospheric density, and the index 
n is weakly dependent upon shielding. However, that method does not work for flux and 
fluence. Instead, we extend this idea by showing that the power law approximation for 
flux J is actually bivariant in energy E as well as density p. The resulting relation is 
J(E,p)~YA (ET )p' n , with A itself a power law in E. This provides another method for 
calculating approximate proton flux and lifetime at any time in the solar cycle. These in 
turn can be used to predict the associated dose and dose rate. 
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1. Introduction 


Studies of space radiation and its effects are concerned with the impact of charged 
species on the functionality and lifetime of human beings as well as scientific 
instrumentation and advanced electronic systems in space. Two aspects of the near-Earth 
space environment are very relevant, particularly in the thermosphere (85km</K500km) 
where Shuttle and International Space Station (ISS) orbits occur. One is the existence of 
energetic proton and electron populations trapped by the Earth’s magnetic field in “Van 
Allen” belts (E.g., Schulz and Lanzerotti, 1974; Spejeldvik and Rothwell, 1985). The 
other is the realization (Jacchia, 1960, 1961) that the properties of the upper atmosphere 
of the Earth are strongly coupled to solar activity, in particular atmospheric density and 
temperature. Throughout the course of the solar cycle, the Earth’s atmospheric neutrals 
expand and contract the thermosphere in response to the behavior of the Sun. Clearly, 
the density in Jacchia’s concept of a dynamic atmosphere couples to the charged-belt 
species as these undergo multiple scattering off the neutrals. That in turn reduces their 
lifetime in the belts (Blanchard and Hess, 1964; Cornwell et al., 1965; Dragt, 1966; Ray, 
1966; Kern, 1994; Pfitzer, 1989; Watts et ah, 1989). 

Therefore, it becomes necessary to understand how atmospheric density per se 
couples to charged-belt population levels as a function of solar activity. This is the 
simplified goal of the present investigation. 

Pfitzer (1989; 1990) has succeeded in developing a reasonable parametric method 
for estimating dose in the thermosphere from atmospheric density. However, the method 
does not work for flux. Inspired by that preliminary investigation, Badhwar and his 
colleagues (1999, 1997, 1996a,b; Golightly et ah, 1996) have examined flight data for a 
correlation between dose and atmospheric density. They have extensively studied and 
analyzed the low-Earth radiation and time lag of the twenty-two year solar modulation of 
the trapped proton radiation exposure inside the Space Shuttle. They have shown that the 
daily trapped-particle dose rate is an approximate power law function of daily 
atmospheric density, thus supporting the Pfitzer model and method. Their further 
analysis of the trapped absorbed dose rate, D, at six fixed locations in the habitable 
volume of the Shuttle exhibits a power law relationship, D = Ap' n , where p is the 
atmospheric density. The index, n, is weakly dependent on the shielding, decreasing as 
the average shielding increases (Badhwar, 1999). 

This present study further examines the AP8 proton flux question and its 
relationship to atmospheric density. It enhances the previous Pfitzer and Badhwar 
density analyses by developing a dynamic trapped-belt proton radiation algorithm that is 
applicable to the ISS and other space flights in the Earth’s thermosphere throughout the 
solar cycle. Although only a very limited range of energies is considered, the method 
addresses several of the shortcomings and over-simplifications in that earlier work. 
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2. Analysis 

The limitations with the original NASA trapped-belt models (Sawyer and Vette, 
1976; Bilitza, 1987) known as AP8 and AE8 have been thoroughly discussed (Watts et 
ah, 1989; Pfitzer, 1989, 1990; Badhwar, 1999). The AP8 model was constructed from 
satellite data in solar cycle 20, a small one compared to more recent events. AP8MIN 
derives from the epoch of 1964, and AP8MAX from that of 1970. The solar radio flux at 
10.7 cm, Fioj, is 150 for AP8MAX and 70 for AP8MIN. These baseline values will be 
adopted here. 

One other promising approach to overcoming the APS model limitations has 
already been produced. It involved the development of a new computer technique known 
as NOAAPRO (Huston and Pfitzer, 1998a, b). This method has since been adapted by 
Singleterry et al. (2004) to enhance the out-of-date AP8 and AE8 models at Shuttle and 
ISS altitudes using the computer program SIREST. 

Since the original AP8 model is readily available elsewhere (Heynderickx et ah, 
2004), it will be used to modify the atmospheric density method of Pfitzer and Badhwar 
by producing a bivariant energy-density algorithm and then compare the result with the 
NOAAPRO-enhanced AP8 model of Singleterry et al. At the outset, AP8 is adopted here 
primarily in order to be consistent with the Pfitzer method. The analysis can be applied 
to other simulation methods such as NOAAPRO and SIREST. Only the omnidirectional 
fluxes are studied in this analysis, noting that the anisotropic nature of these has been 
discussed by Watts et al. (1989). 

Upon examining the proton flux data from the AP8 model program, and in view 
of the overall problem as studied for more than 40 years (Pfitzer, 1989, 1990), several 
observations can be summarized. 

1 . The atmospheric density decreases rapidly as the altitude increases. 

2. Both integrated and differential flux for protons increases as the altitude 
increases (or density decreases) for both solar maximum and minimum cycles. 

3. The fluxes in general are smaller at solar max than at solar min, at least for 
low altitudes and low energies. 

4. However, the difference in the proton flux is wider at low energy than for the 
higher-energy protons at the same altitude. 

5. The flux decreases as the energy increases. 

These observations prompt the idea that the proton flux / can be expressed, 
empirically at least, as a function of two variables, density (altitude) and energy, namely 
J(E,p) ovJ(E,h). 


Solar Modulation of Atmospheric Density 

The altitude and density relationship has a long history (Jacchia, 1960, 1961; 
Harris and Priester, 1962, 1963). The form to be used here has been described recently 



by several authors (Pfitzer, 1989, 1990; Watts et al., 1989). A simple parameterization of 
the US Air Force Model made by Pfitzer (1989, 1990) is 

p = p 0 exp{-(h-120) / [M(h-103) I/2 ]}, (1) 

where the solar-cycle modulation term M is expressed as 


M= 0.99 + 0.518 { (F io .7 + F bar ) / 110 } ,/2 , (2a) 

F = F io.7 + Fbar ■ (2b) 

In (1) p denotes the atmospheric density, p 0 is assumed to be p Q - 2.7 x 10' 11 g/cm 3 , h is 
the altitude in km, Fw.? is the instantaneous value of the solar radio flux at 10.7 cm, and 
Fbar is the weighted value of F 10.7 for averaging, such as three prior solar rotations. 

The density in (1) is a multi -variant function of h and F,oj. Similarly, the AP8 
proton flux J is a multi-variant function of h and energy E. The problem at hand, then, is 
to generate the multi-dimensional surface of J as a function of E and h or p. By selecting 
an altitude h and emulating the solar cycle with a solar radio flux F )o. 7 , the modulated 
proton flux J follows as a function of energy E. Dynamically speaking, furthermore, all 
of these variables are functions of time t. 

A “carpet” plot (in the sense of Pfitzer, 1989, 1990) is merely a projection of 
these surfaces onto a two-dimensional graph. This can be obtained for the solar-cycle 
terms altitude h -f(p) and flux F 10.7 = g(p) as a function of atmospheric density pin (1) 
and (2), by taking the inverse of (1) for constant surfaces of F 10.7 and h respectively. The 
result is provided in Fig. 1. As pointed out by Pfitzer, the trapped protons have a very 
slow response time to dynamic changes in atmospheric density p(t). Therefore, the 
values of F 10.7 and Fbar are assumed identical, whereby F = F , 0.7 + Fbar = 2 F 10.7 in (2b). 
As stated earlier, the values of F,o. 7 used for solar max and solar min in the following 
regression analysis become 150 and 70 respectively. 

Regression Algorithm 


Several methods and approaches are available to generate a semi-empirical 
formula for proton flux J as a bivariant function of density p (altitude h ) and energy E. 
The method adopted here is taken from theoretical nuclear physics (Lodhi and Waak, 
1975, 1976) based upon a polynomial regression analysis. It is used to determine the 
functional relationship between fluxes at solar maximum and minimum. Since the proton 
lifetime x is determined by energy losses primarily due to multiple Coulomb scattering 
from charged constituents in the upper thermosphere, as well as some neutral scattering, 
it is a function of time t, proton energy E, and atmospheric density p or altitude h. That 
is, x = x(E, p l , t) or x = x(E, h, t) since the atmosphere expands and contracts at different 
times during the solar cycle. 
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Utilizing the regression analysis technique, one keeps the regression coefficient 
greater than 90%. A ratio J max /J m i n or J m i„/J max is generated for differential proton fluxes 
at solar max and solar min as a function of density (and 1/density) for energies between 
30 and 300 MeV from AP8 model data within 300 to 600 km. Next this ratio is fitted to 
some polynomial ranging from linear to fifth-order in p and 1/p for proton energies of 30, 
35, 40, 45, 50, 60, 70, 80, 90, 100, 125, 150, 175, 200, and 300 MeV. Note that this ratio 
is non-linear in p or 1/p. One also finds that polynomials of higher-order than second 
result in a better fit for given energy than the second-order. However, a common 
expression for the entire energy range gets worse than that of the second-order. This 
observation limits the method to be confined to a second-order expression for the flux 
ratio as a function of p max or p l max for all energies. The regression analysis works well 
within the energy and altitude range adopted. For other ranges of approximation one has 
to be careful and do the analysis again, piecewise fitting the entire dynamical range. 

Density Dependence 

First the flux ratio is generated as a function of p in the quadratic form. The 
coefficients of p vary while progressing from one energy to the next. The flux ratio can 
then be written in the following fashion: 

d min/ J max ~ dp 3 " bp 3 ” C , ( 3 ) 

where a, b, and c are proton energy coefficients. Naively, one might hope that these are 
constant coefficients. However, one discovers that a, b, and c are functions of energy E 
for the 15 different energy values chosen. 

The next step is to find a common expression for this ratio for all energy values. 
That is achieved by obtaining a functional relation for the coefficients in (3) as a function 
of energy, again by the method of regression. One obtains four relations for p ranging 
from linear to fourth-order in energy. One also finds that expressions for coefficients a, 
b, and c cannot be of the same polynomial order in reproducing the flux ratio. The best 
fits found for two different energy ranges are as follows: 

J min/ J max = faE 2 + d]E + do)l e p 2 + faE* + b 3 E 3 + b 2 E 2 + b}E + bo)leP 

+ (cjE 4 + cjE 3 + C 2 E 2 + cjE + co)ie 1 30 < E < 60 MeV (4) 


Jmin/Jmax (djE + do) he P 3 " (b jE + bo) he P 


+ (c 2 E/ + ciE + c 0 )he ■ 70 < E< 300 MeV ( 5 ) 
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The coefficients within parentheses are different for le (low energy) and he (high energy). 
The actual coefficients and a check for the accuracy are given below under Results. 

Inverse Density Dependence 

In contrast to (3), the inverse algorithm can be derived. It is known that the 
atmospheric densities decrease as the altitude increases or the reciprocal of the density 
increases as the altitude increases. The flux variation with respect to the inverse of the 
density must convey a direct relation to the variation of the altitude. One must therefore 
search for a similar expression giving the flux ratio as a function of the inverse density. 
After several trials the best-fitted function thus obtained is given in the form: 

Jmcu/Jmin ~ (&4 E? + ct 3 E 3 + 1 12 E 2 + ctjE + cio)p~ Jr (P 3 E 3 + b 2 E 2 + bjE 

+ bo)p I + ( C 2 E 2 + c,E + Co) (6) 

for all energies E under consideration. This expression is further approximated by 
writing the coefficients in the exponential form, thus yielding: 

JmJJmin = Ae^p 2 + BeP E p 1 + Ce rE (7) 

for all proton energy ranges from 30 to 300 MeV. 

Results at Solar Max 

The two algorithms (3) and (7) are now compared at solar maximum. The 
resultant semi-empirical formula for the flux ratio in (3) as a function of p (in units of 10 " 
15 g/cm 3 ) is given by: 


(Jmin/JmJle = (-3x1 O^E 2 + 1x1 O' 5 E - 8.0x1 )p 2 

+ (7xl0' I0 E? - 5x 10' 7 E 3 + lxl O' 4 E 2 - 1.4xlO' 2 E + 0.695) p 
+ (-2x10'“ E 4 - 5x1 O' 8 E 3 + 4x10' 5 E 2 - 8.6x1 O' 3 E + 1.897) , 

30 <E <60 MeV ( 8 ) 

(Jmin/Jmadhe = (2x1 O' 6 E + lxl0' 4 )p 2 + (-7x1 O' 4 E + 0.181)p 
+ (6x1 O' 6 E 2 - 3. 7x1 O' 3 E + 1. 66) . 

70 <E <300 MeV (9) 
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To check formulas (8) and (9) for 7 W( „, an example at 400 km and proton energy 
of 100 MeV is taken, and the results summarized in Table 1. Let us define 

f(p) = ap 2 + bp + c (10) 

on the right-hand-side of (3), (8), and (9). 

At altitude 400 km (Heynderickx et al., 2004), the AP8 model in SPENVIS gives 

Pmax = 3.8xl0' l5 g/cm 3 

Pmin - 9.57x1 0' 1 6 g/cm 3 

JnJ }00MeV) = 2. 79x1 a 2 cm 2 s’ MeV 1 (SPENVIS) 


For this density at solar max one obtains f(p max ) =f(3.8) in (10) and J m i„/J n tax = 1 ■ 78 in 
(3). It then follows from (3) that 

J,J Algonthm) = J m J SPENVIS) f(p max ) = (2.79x1 O' 2 cm’ 2 s' 1 Me V 1 ) (1.78) 

= 4.9 5 5 4x1 O' 2 cm 2 s' MeV 1 . (1 1) 

On the other hand, AP8 (Heynderickx et al., 2004) gives 

Jmin (SPENVIS) = 5. 15 1x1 O’ 2 cm 2 s‘ MeV 1 . (12) 

Comparison of (1 1) and (12) yields a difference of 0.1 83x1 O’ 2 cm' 2 s' 1 MeV 1 with an error 
of 3.5%. These are summarized in Table 1 . 

Next, the resultant semi-empirical formula for the flux ratio JmcJJmm in (7) as a 
function of 1/p (in units of 10 +I5 cm 3 /g ) is determined by regression analysis to have the 
form: 

Jmax/Jmin = -0.0241 e° 0007E p 2 + 0. 1 966e° 0007E p 1 

+ O.3208e +00032E (13) 

for all energies between 30 and 300 MeV. 

Following the same procedure used in (10) and (1 1), we can define 


g(p l ) = Ae aE p 2 + Be /3E p‘ + Ce yE 


(14) 



for the right-hand-side of (7), (8), (9), and (13). One determines that g(p 1 max) = 0.4912 in 
(14). The resultant J min , the difference from SPENVIS, and the error are summarized in 
Table 1. 

For further comparison, from expression (13) for J max {cm‘ s' 1 MeV ! ) a differential 
flux is calculated and contrasted with the AP8 data in Fig. 2, derived from SPENVIS 
(Heynderickx et al., 2004) for an ISS orbit of 350-478 km altitude and inclination 51.65°. 
Fig. 2 is a two-dimensional projection of the three-dimensional surface J(h,E) at various 
selected altitude h. The NOAAPRO results in SIREST are shown in Fig. 3 and Fig. 4 
along with the algorithm at 400km and 500km altitudes, for solar max with Fw .7 = 150 in 
the algorithm. 

Results Midway between Solar Max and Solar Min 

Method 1. By varying the solar-cycle modulation parameter Min (2), one obtains 
a different atmospheric density model in (1). This can be accomplished by changing F 10.7 
and Fbar whereby a different value of density p is obtained, either from (1) or from the 
carpet plot in Fig. 1. Upon insertion of the new value of density p, a proton differential 
flux follows from (3) and (7). The baseline regression algorithms (3) and (7) assume 
F10.7 - 70 and F10.7 = 1 50 in Fig. 3 and Fig. 4 for solar min and max respectively. In 
order to determine the proton differential flux mid-way through this same adopted cycle, 
the solar flux becomes F )o. 7 ~ 110 whereby F = 220 in (2b). The resulting proton 
differential spectrum is shown in Fig. 5 for 400km in Fig. 6 for 500km. 

Method 2. Given a proton flux J at either solar maximum or minimum, such as 
the algorithm in (3) and (7), then an interim flux is determinable as a linear time- 
interpolation, 

J(E,h,T)~J max (E,h,r)(l-r) + F J min (E,h,z) , (15) 

or alternatively, 

J(E,h,T)~J min (E,h,T)(l-r) + r J max (E,h,r) , (16) 

where r(E,p, x) is the dimensionless fraction 


F (E, p, x ) = r ; Tmax(E) . (17) 

*min(E) - ?max(E) 

The lifetime x is assumed to be limited to one solar cycle or 1 1 years. 

To calculate the desired intermediate proton flux J(E,h,r) at a time between solar 
maximum and solar minimum using Method 2, the right-hand-side of (17) represents the 
interpolation fraction r of the solar cycle since last solar minimum. Then either of (1 5) 
and (16) gives the interim flux in this approximation. Substituting (3) and (7) into (15) 
and (16) respectively, one has 
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J(E, h,r)~ J max (E, hr) \(1 -JT) + r ( ap +bp + c ) 


(18) 


J(E,h,T)~J min (E,h,T)p-r) + r(Ae aE p 2 +Be^ E p l +Ce yE ) J . 


(19) 


The various coefficients in (18) and (19) are given in (4)-(5) and (8)-(9). Further study is 
planned to conduct an error analysis between Method 1 and Method 2. 


3. Conclusions 

The proton differential flux has been expressed, empirically, as a bivariant 
function of density (altitude) and energy, broken into two ranges of proton energy, viz., 

30 to 60 MeV and 70 to 300 MeV. The corresponding expression in terms of inverse 
density is relatively compact and works for the entire range of proton energy, 30 to 300 
MeV. From this analysis it is observed that the proton differential flux has a nonlinear 
dependence on energy and density (altitude). The flux ratio has been expressed in a 
semi-empirical form given by (3) and (7). It has been compared with AP8 model data as 
shown in Fig. 2 for Shuttle and ISS altitudes of current interest. An additional 
comparison with NOAAPRO is given in Fig. 3 and Fig. 4. An illustration of the 
algorithm mid-way through the adopted solar cycle is produced in Fig. 5 and Fig. 6. 
Finally, the algorithm provides a simple means for calculating the trapped-belt proton 
flux when the Fwj solar modulation flux is 200 or greater. The analysis thus avails itself 
to other more prominent solar sycles where AP8 is not valid. However, a thorough error 
analysis will be required in the future in order to determine the general limitations of this 
regression-analysis algorithm. As a concluding remark, the selection of solar flux Fjo.7 is 
a matter of convention due to its known correlation with sunspot number. The physics of 
the thermosphere is not completely understood and there is current interest in the extreme 
ultraviolet parameter Ejoj. Should another modulation factor be found, the regression 
analysis presented here can be modified to accommodate it. 
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Figure 1 Caption 

Carpet plot of the solar-cycle terms altitude h and modulation flux F10.7, as a function of 
atmospheric density in Equations (1) and (2). 

Figure 2 Caption 

Graph of proton differential flux versus energy at various Shuttle and International Space 
Station altitudes, comparing the present algorithm with predictions of AP8 at solar 
maximum. 
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Figure 3 Caption 

Graph of proton differential flux versus energy at 400 km altitude. Proton flux models 
AP8MAX, SIREST/NOAAPRO, and the algorithm presented here are compared at solar 
maximum with F 10.7 = 150. AP8MIN is also given. 

Figure 4 Caption 

Graph of proton differential flux versus energy, like Figure 3, except at 500 km altitude. 
Proton flux models SIREST/NOAAPRO, and the algorithm presented here are compared 
at solar maximum with F 10.7 - 150. AP8MAX is not illustrated since it is essentially 
identical to SIREST at this altitude. AP8M3N is also given. 

Figure 5 Caption 

The same graph as in Figure 3, except with the algorithm initialized for half-way through 
the assumed solar cycle assuming F I0 ,7 - F bar = 1 10 in Eq. (1) and (2). 

Figure 6 Caption 

The same graph as in Figure 4, except with the algorithm initialized for half-way through 
the assumed solar cycle assuming F !0 .7 = F bar ~ 1 10 in Eq. (1) and (2). 
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Graph of proton differential fbx vs. energy at various Shuttle and International Space Station altitudes 
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Graph of proton differential flux vs. energy at 400km 
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Figure 3 
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Figure 4 
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Energy (MeK/) 


Figure 5 
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Figure 6 
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Table and Table Caption 


Table 1 

Comparison of J min (Aigor “ hm) with J mi „ (SPENVIS) for proton energy of 100 MeV at 400 km. 


SPENVIS Parameter 

Value 

Value 

Pmax 

3.8x1 O' 15 


Pmin 

9. 5 7x1 O' 16 


r (SPENVIS) 
Jmax 

2.79x1 O' 2 


t (SPENVIS) 
J min 

5.1 5x1 O' 2 


r (SPENVIS) / r (SPENVIS) 
Jmin ' Jmax 

1.78 


Comparison with Algorithm 

Algorithm (3) 

Algorithm (7) 

r (Algorithm) 
Jmin 

4.955 

5.680 

Difference from J min (SPENVIS > 

0.183 

0.53 

% Error 

3.5 

10.3 


Units of p are in g/cm 3 and those of J are cm 2 s’ 1 MeV' 1 . 
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